// This is core/vnl/algo/tests/test_functions.cxx
#include <iostream>
#include <iomanip>
#include <cmath>
#include <testlib/testlib_test.h>
//:
// \file
//
// Test function results for chi-squared cumulative density function.
// The values are those calculated from this function on solaris, and
// agree (to 3sf) with those from a mathematical table.
//
// Each row is for a different dof, from 1 through 96 in increments of
// 5 (20 rows).
//
// Each column is for a different chi-squared, from 0 though 180 in
// increments of 20 (10 columns).
//
// For readability reasons, the rows are split in 2x 5 rows, with an extra
// indentation for the second half of each row.
#include <vnl/algo/vnl_chi_squared.h>
double cdf_baseline[] =
{
  0.0, 0.9999922488859869, 0.9999999997449120, 0.9999999999999905, 1.0,
     1.0, 1.0, 1.0, 1.0, 1.0,
  0.0, 0.9972306042844884, 0.9999995444850495, 0.9999999999549898, 0.9999999999999964,
     1.0, 1.0, 1.0, 1.0, 1.0,
  0.0, 0.9546593186683575, 0.9999642248743472, 0.9999999907278384, 0.9999999999985243,
     0.9999999999999998, 1.0, 1.0, 1.0, 1.0,
  0.0, 0.7797793533983011, 0.9992214099174926, 0.9999994766265833, 0.9999999998335990,
     0.9999999999999654, 1.0, 1.0, 1.0, 1.0,
  0.0, 0.4787387426184652, 0.9925632227010702, 0.9999872283117348, 0.9999999919251898,
     0.9999999999971140, 0.9999999999999993, 1.0, 1.0, 1.0,
  0.0, 0.2084435236051257, 0.9609880071450072, 0.9998323023504486, 0.9999997903577246,
     0.9999999998716507, 0.9999999999999505, 1.0, 1.0, 1.0,
  0.0, 0.0641963723612628, 0.8709596508694845, 0.9986502140214748, 0.9999966874714471,
     0.9999999965356298, 0.9999999999979206, 0.9999999999999991, 1.0, 1.0,
  0.0, 0.0142776135970496, 0.7029716020753263, 0.9927297837948102, 0.9999652953758086,
     0.9999999382046935, 0.9999999999423534, 0.9999999999999658, 1.0, 1.0,
  0.0, 0.0023551011658856, 0.4850483796984979, 0.9720590475678458, 0.9997436372936986,
     0.9999992259610407, 0.9999999988798406, 0.9999999999990349, 0.9999999999999994, 1.0,
  0.0, 0.0002957368080061, 0.2793886568739744, 0.9194309789050203, 0.9986010626587202,
     0.9999928735024521, 0.9999999840351245, 0.9999999999800138, 0.9999999999999838, 1.0,
  0.0, 0.0000289361257476, 0.1331225579748673, 0.8182400707526587, 0.9941486753062340,
     0.9999500186862800, 0.9999998271108139, 0.9999999996858998, 0.9999999999996467, 0.9999999999999998,
  0.0, 0.0000022535340508, 0.0524807132282665, 0.6671309159544766, 0.9806610209037141,
     0.9997251552759223, 0.9999985364764852, 0.9999999961469799, 0.9999999999940054, 0.9999999999999938,
  0.0, 0.0000001352148233, 0.0172069984403675, 0.4878147556882595, 0.9481646213453065,
     0.9987862811130113, 0.9999900872908165, 0.9999999622464486, 0.9999999999188310, 0.9999999999998879,
  0.0, 0.0000000073722077, 0.0047274255385307, 0.3154587502878544, 0.8846964150660444,
     0.9956070677677689, 0.9999452153818948, 0.9999996987525844, 0.999999999105923, 0.9999999999983556,
  0.0, 0.0, 0.0010974342908522, 0.1790828046260967, 0.7825800096649158,
     0.9867347062121498, 0.9997487615314837, 0.9999980101056154, 0.9999999918570519, 0.9999999999800666,
  0.0, 0.0000000000116404, 0.0002170786956947, 0.0890129923177431, 0.6453492639226499,
     0.9660451059031501, 0.9990298866149123, 0.9999889634183597, 0.9999999378147725, 0.9999999997975395,
  0.0, 0.0, 0.0000368857904673, 0.0387583535447800, 0.4894571445703249,
     0.9252463664134718, 0.9968048573119753, 0.9999479616312844, 0.9999995969280049, 0.9999999982561550,
  0.0, 0.0000000000000097, 0.0000054252284960, 0.0148195154453572, 0.3381826989258588,
     0.8564977678812210, 0.9909185509763581, 0.9997891030141784, 0.9999977585673097, 0.9999999871266845,
  0.0, 0.0, 0.0000006956282422, 0.0049930250351178, 0.2116235676678424,
     0.7565873548899793, 0.9774854404134016, 0.9992580809522048, 0.9999892041010620, 0.9999999177818010,
  0.0, 0.0000000000000002, 0.0000000782693028, 0.0014883025723802, 0.1195827092276646,
     0.6303318279959211, 0.9508233262755118, 0.9977140253664275, 0.9999545743054137, 0.9999995418724694
};

int test_functions()
{
  int idx = 0;
  for (unsigned int n=1; n<100; n+=5)
  {
    for (unsigned int chisq_int = 0; chisq_int < 200; chisq_int += 20)
    {
      const auto chisq = static_cast<double>(chisq_int);
      const double cdf = vnl_chi_squared_cumulative(chisq,n);
      const double err = std::fabs(cdf - cdf_baseline[idx++]);
      std::cout << "vnl_chi_squared_cumulative(" << chisq << ',' << n << ')';
      TEST_NEAR(" CDF", err, 0.0, 2e-15);
      if (err >= 2e-15)
        std::cout << "Error = " << std::setprecision(16) << err << std::endl;
    }
  }

  std::cout << "cdf(7.88,1) = " << vnl_chi_squared_cumulative(7.88,1)
           << " should be about 0.995\n"
           << "cdf(14.8,12) = " << vnl_chi_squared_cumulative(14.8,12)
           << " should be about 0.75\n"
           << "cdf(10.1,19) = " << vnl_chi_squared_cumulative(10.1,19)
           << " should be about 0.05\n"
           << "cdf(39.3,40) = " << vnl_chi_squared_cumulative(39.3,40)
           << " should be about 0.50\n"
           << "cdf(109.1,100) = " << vnl_chi_squared_cumulative(109.1,100)
           << " should be about 0.75\n";

  // rand() is not always a good random number generator,
  // so use the following congruential random number generator - PVr
  static auto sample_seed = (unsigned long)std::time(nullptr);

  double hist1[20];
  for (double & i : hist1)
  {
    sample_seed = (sample_seed*16807)%2147483647L;
    double u = double(sample_seed)/0x7fffffff; // 0x7fffffff == 2147483711L
    i = 10.0+20.0*(u-0.5); // uniform in the interval 0 - 20
  }
  double chisq = 0;
  for (int i=0; i<20; i++)
  {
    std::cout << i << ' ' << hist1[i] << std::endl;
    double delta = hist1[i] - 10.0;
    chisq += delta*delta/(hist1[i] + 10.0);
  }
  std::cout << "cdf(" << chisq << ",20) = "
           << vnl_chi_squared_cumulative(chisq,20)
           << " so P(same dist) = " << (1.0 - vnl_chi_squared_cumulative(chisq,20))
           << std::endl;
  return 0;
}

TESTMAIN(test_functions);
